Understanding and controlling A/^-dimensional quantum walks via dispersion relations. 
Application to the 2D Grover walk: Diabolical points and more 
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Abstract 

The discrete quantum walk in N dimensions is analyzed from the perspective of its dispersion relations. This allows under- 
standing known behaviors as well as designing new ones when spatially extended initial conditions are considered. This is done 
by deriving wave equations in the continuum, which are generically of the Schrodinger type, and allow devising interesting 
behaviors, such as the ballistic propagation without deformation or the generation of almost flat probability distributions, 
what is corroborated numerically. There are however special points where the energy surfaces display conical intersections (the 
so-called diabolical points) and, near them, the dynamics is entirely different, similar to that of massless particles or of electrons 
in graphene. Applications to the two-dimensional Grover walk are presented. 

PACS numbers: 03.67.-a, 42.30.Kq 
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I. INTRODUCTION 

Quantum walks (QWs) are becoming widespread in 
physics. Originally introduced as quantum versions of 
random classical processes [1] and quantum cellular au- 
tomata [2], they have been deeply investigated, specially 
in connection with quantum information science [3|49|, 
and have been recently shown to constitute a univer- 
sal model for quantum computation p^, lUJ]. Remark- 
ably, along the last decade this simple quantum diffusion 
model has found connections with very diverse phenom- 
ena such as optical diffraction and other wave phenom- 
ena lH, 13 1, quant uni games [l^ . Anderson localization 
and quantum chaos jl5l4l8| . and has even been identi- 
fied as a generalized quantum measurement device [19j]. 
Moreover, through proposals for and actual implemen- 
tations QWs have found connections with the dynamics 
of quite different physical systems such as linear optical 
systems and optical cavities [20- 28j , optical lattices [29- 
[32I . trapped ions 33l-[35| , and Bose-Einstein condensates 
[36] among others. Along these recent years the exper- 
imental capabilities in performing QWs have increased 
steadily and nowadays there are devices that, e.g., im- 
plement two dimensional QWs with a large number of 
propagation steps [37.]. For recent reviews on the field 
see (sil. [sot- 

Recently, some of us have been insisting in the interest 
of analyzing the dynamics of coined QWs from the per- 
spective of their dispersion relation [l3|, lioj, w hich has 
been also used for purposes different to ours [4l|, |42| . We 
think this is an interesting approach that is particularly 
suitable for non localized but extended initial states, i.e., 
when the initial state is described by a probability distri- 
bution extending on a finite region in the lattice. When 



this region is not too small, the evolution of the QW can 
be derived to a good degree of approximation from famil- 
iar linear differential wave equations. Certainly the use 
of extended initial conditions has been rare up to now 
[4^ , but we think it is worth insisting in that they allow 
to reach final probability distributions which, interest- 
ingly, can be tailored to some extent. It is also worth 
mentioning that extended initial conditions also make a 
connection with multiparticle quantum walks for nonin- 
teracting particles [43 1. 



In [40] we have applied the above viewpoint in studying 
alternate QWs [45] in N dimensions. Alternate QWs are 
a simpler version of multidimensional QWs than the stan- 
dard one as they only require a single qubit as coin 
Q instead of using a 2A^-qudit as coin that is used once 
per time step, as originally defined in [46] (the price paid 
for this simplification of coin and coin operators is that 
the coin operator must be applied N times per time step, 
alternating with the displacement operator which must 
be applied along each one of the N dimensions). Here 
we shall carry out, from the dispersion relation view- 
point, the study of multidimensional QWs in their stan- 
dard form. Below we make a general treatment valid for 
N dimensions and illustrate our results with the special 
case of the the Grover two-dimensional QW, which has 
received most of the theoretical attention paid to multidi- 
mensional QWs [46145 7( . We shall derive continuous wave 
equations, whose form depend on the probability distri- 
bution at time t = 0, which allow to choose a predeter- 
mined behavior for the time distribution. We also show 
that the A^-dimensional QW may present special points 
in their dispersion relation called diabolical points [58- 



60] that determine a dynamics similar to that of massless 
Dirac fermions or of electrons in graphene [6ll-[65|. a re- 
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suit we already found in j4o[ and that we treat here with 
mathematical detail. We want to insist in the universal- 
ity of the diabolical points, which are conical intersec- 
tions in the dispersion relation, as they appear in physics 
in very different systems: Apart from the already men- 
tioned graphene electrons they happen in the spectra of 
polyatomic molecules [gg], quantum triangular bil- 
Hards conical refraction in crystal optics [Hoj . 
optical lattices [sH or acoustic surface waves [69|, just to 
mention a few. In all cases the remarkable physical prop- 
erties exhibited by all these systems is directly associated 
with the existence of the diabolical points. We will see 
that the Grover QW in 2D also exhibits this remarkable 
feature, as Alternate QWs do indeed [40], but at differ- 
ence from all the above mentioned systems, here we are 
dealing with discrete dynamical systems. 

The rest of the paper is organized as follows. In the 
next section we introduce the TV- dimensional QW and 
solve it formally. Then in the third section dispersion 
relations and diabolical points are presented for the spe- 
cial case of the two-dimensional Grover QW. In section 
IV general continuous wave equations are derived for the 
case of spatially extended initial conditions, and then the 
2D Grover walk is specifically treated. In the last section 
some conclusions and outlooks are presented. 

II. N-DIMENSIONAL DISCRETE QUANTUM 
WALKS. GENERALITIES 

In the A/'-dimensional QW the walker moves at discrete 
time steps t G N across an A/'-dimensional lattice of sites 
X = (xi,...,XAr) G Z^. The walker is endowed with 
a 2A/'-dimensional coin which, after a convenient unitary 
transformation, determines the direction of displacement. 
The Hilbert space of the whole system (walker+coin) has 
then the form 

n = np^nc. (i) 

where the position space, Hp, is spanned by 
{|x) = . . . ,XAr) : G Z; a = 1, . . . , A/"} ((x |x') = 
and the coin space, He, is spanned by 2N or- 
thonormal quantum states {\ar}) : a = 1, . . . , A/"; = ±}. 
Note that a is associated with the axis and r] with the di- 
rection. For example, in the popular ID QW (A^ = 1) we 
would have just |1+) and which are the equivalent 
to the |R) and |L) (for right and left) states commonly 
used in the literature. We introduce this apparently com- 
plicated notation in order to deal easily with an arbitrary 
number of dimensions. 

The state of the total system at time t is represented 
by the ket \^l^t)^ which can be expressed in the form 

N 

i^*) = EEE<'' w®K)' (2) 

X a—1 77=± 



where the projections 

= (KI®(x|)|f*), (3) 

are wave functions on the lattice. We find it convenient 
to define, at each point x, the following ket 

N 

l^x.*) = (x|v*) = 5]EC.f K), (4) 

which is an (unnormalized) coin state, so that ?/^^'^ = 

{(^T] I V^x,t). As 1^ = \{{ar^\ (x|) \iljt)f is the prob- 
ability of finding the walker at (x, t), and the coin in 
state the probability of finding the walker at (x, t) 
irrespectively of the coin state is, then, 

^-.*=EEK''r=(^x,ti^x,t), (5) 

a=l 77=± 

where we used the fact that X]^=i ^n=± 

identity in He- Clearly ^^Px,t = 1 because |x) (x| 

is the identity in Hp. 

The dynamical evolution of the system is ruled by 

\^t+i) = U\4't), (6) 
where the unitary operator 

U = Do(^i^C^ (7) 

is given in terms of the identity operator in Hp, /, and 
two more unitary operators. On the one hand C is the 
so-called coin operator (an operator in He), which can 
be written in its more general form as 

N 

c= J2 E ca'K)K'|, (8) 

a,a' = l rj,r]'=^ 

where the matrix elements C^]^^, = C |q^J^/) can be 
arranged as a 2A^ x 2N unitary square matrix C. On the 
other hand D is the conditional displacement operator in 
H 

AT 

^ = E E E 1^ + ^^-) (^1 ^ 1^^) Kl . (9) 

X a=l 77=± 

where Uq, is the unit vector along direction Xq,; note that, 
depending on the coin state lofr^), the walker moves one 
site to the positive or negative direction of if 77 = + 
or 77 = — , respectively. 

Projecting (j6j) onto (x| and using (j3j) and ([7j)-(j9j) we 
get straightforwardly 

N 

l^x,t+i> = E E I"-?) Kl c |^x-^u„,t) , (10) 
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which further projected onto (a^yl leads to 

N 

:,t+l ■ ^ a'r}' TX—riUa.t' 



a' = l 77'=it 



Equation (p!Q|) . or equivalently (pTj) . is the NDQW map 
in position representation; it shows that, at each (dis- 
crete) time, the wavefunctions at each point are coher- 
ent hnear superpositions of wavefunctions at neighboring 
points at previous time, the weights of the superposition 
being given by the coin operator matrix elements C^J^^, . 
Given the linearity of the map and the fact that it is 
space-invariant (C^;^/ do not depend on space) a useful 
technique here is the spatial Discrete Fourier Transform 
(DFT), which has been used many times in QW studies. 
First we define the DFT pair 



^/;k,i)^Ee""'"i^x,t), 



1^: 



_ikx 



(27r) 



(12) 
(13) 



where k = (/ci, . . . , /cat) and ka G [— 7r,7r] is the (quasi- 
) momentum vector [70]. According to (j4]) we get 



N 



E E I".) , 

a=l 77=± 



(14) 
(15) 



Applying the previous definitions to the map (pTj) . and 



using Y^^e- 
get 



-ikx 



1^: 



X-77U«,t/ 



-ir]'k-Ua 



'0k, t / we readily 



(16) 



where we defined a coin operator in the quasi-momentum 
space 

N 

C'k^EEK)KIC'e-''"", (17) 
ka = k • Uq,, whose matrix elements read 

Ki Ck K) ^ {Cur^i, = e-^^'-c^'- (18) 

Projection of (p!6|) onto and using (|17ll8p leads to 

N 

fei = E E ^-^'''■""C':;'^)'''- (19) 



Hence the nonlocal maps (|10|lip become local in the mo- 
mentum representation (|16ll9p . This allows solving for- 
mally the QW dynamics very easily because map ([16] 
implies 



'0: 



(20) 



and hence the eigensystem of (7k (or of Ck in matrix 
form) is most useful for solving the problem, as we do 
next. 

As the operator Ck is unitary (it must be so according 



to ([20]) ) its eigenvalues |Aj^*^ | have all the form 



exp 



(-i^L-') 



with 



(s) 

real. Apart from the eigenvalues 
we will need to know its eigenstates | |0k / 1 

Once the eigensystem of Ck is known, implementing 
([2Q|) is trivial: Given an initial distribution of the walker 
in position representation |'0x,o) we compute its DFT 

'^k,o) via ([12]), as well as the projections 



so that 
get 



V^k,0 



/k 



is) 
k 



k,0 



(21) 



Now recalling ([20]) we 



'k,t 



2N 

E- 



(22) 



where we used A^*^ = exp ^— icj^*^^ , while in position 
representation we get, using ([T3|) . 



I^x,*) = 



2N 

E 



I 



(2.)"' 



i ( k-x— a; 



^/k 



(23) 
(24) 



Hence the QW is formally solved: all we need is to com- 
pute the eigensystem of Ck and the initial state in recip- 
rocal space ?/ik,o ) 5 which determines the weight functions 



fl'^ through 



Equation ([23]) shows that the QW dynamics corre- 
sponds to the superposition of 2A^ independent walks, 
labeled by s. According to (|24|) the co'^^^'s are the fre- 
quencies of the map, each of which defines a disper- 
sion relation in the system {2N in total). Note as 
well that what we have done in the end is to decom- 
pose the QW dynamics in terms of plane waves. In 
particular, if f^^ = J^^^ (k — ko), what means that 



?/^k,o / is different from zero only for k = ko. 



[ZTT) exp 



i (ko 



X 



which is an unnor- 



malizable plane wave and thus unphysical. 

In order to avoid possible confusions, to conclude this 
initial part we state that the ordering of the coin base 
elements we will be using in the matrix representations 
of operators and kets is |1+) , |1_) , . . . |A^+) , 
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III. DIAGONALIZATION OF THE 2D GROVER 
MAP: DISPERSION RELATIONS AND DIABOLI- 
CAL POINTS 

In this Section we illustrate the general treatment in- 
troduced above in the special case of the 2D Grover walk. 
This is a special case of the Grover walk in N dimen- 
sions, which is defined by a matrix coin C (the so-called 
Grover coin of dimension 2A^) of matrix elements given 
by C^;^/ = 1/N — 5a,a'^r],r)' • In the present case {N = 2), 



^=2 



/-II 1 1\ 
1 
1 



-1 1 1 
1 -1 1 

1 1 -1/ 



(25) 



The corresponding matrix Ck ([18]) , with k = (/ci, /C2), has 
the form 





/- 




g-i/Ci 




g-i/ci\ 


1 










2 




g-i/C2 




_g-i/C2 


g-i/C2 




V 


gi/C2 


gi/C2 


gi/C2 


_gi/c2 y 



(26) 



whose diagonalization yields the eigenvalues Aj^ 
exp ^— i^k^^^ with 



is) 



,(1) 



(2) 

: TT + r^k, ^k 



O (3) 



(4) 



0, (27) 




1.D.0 



Figure 1: (Main figure): Dispersion relation sketch around the 

diabolic point. The value of cj/tt is plotted for the different 

sheets (i.e., different eigenvalues). The middle (fiat) surface 

corresponds to uj^^^ = tt, the surface with values ranging from 

1 to 2 corresponds to ojI}\ The bottom plane is associated 

with (jj^^ and the surface with values ranging from to 1 

corresponds to c<;^^\ The smaller inset shows a detail of cj^^'* 
(2) 

and cj^ around the origin. 



where 



- (cos ki + cos /C2) 



G [0,^]. 



(28) 



Note that adding a multiple integer of 27r to any of the 
cj's does not change anything because time is discrete 
and runs in steps of 1, see (|24|) . 

We see, according to (|27|) . that the last two eigenval- 
ues A^^'^^ = exp i^k"^'^^^ = —1,1 do not depend on k, 

and the first two eigenvalues A^^'^^ = exp ^— i^k^'^^^ = 
— exp (^il^k) are complex-conjugate of each other, hence 
we can choose l^k ^ [0, tt] without loss of generality. A 
plot of the dispersion relations (|27|) is given in Fig. 1, 
where five special points (the origin and the four cor- 
ners, see also Fig. [2]) are observed: there we see 3-fold 
degeneracies. Moreover the frequencies il^k have a con- 
ical form (diabolo-like) at those degeneracies and this is 
why we call them diabolical points, in analogy to other 
diabolos found in different systems (58I462L I66l469j . At 
the central diabolical point the degeneracy is between 



UJ- 



(s=l,2,3) 



, while at the corners it is between u-^ 



(s=l,2,4) 



For the sake of later use we note that the frequency l^k 
reads, close to the diabolical point k = 0, 



71 



cos (2(9) 
48^2 



0{k^ 



(29) 



where k = k {cos 6 , sin 6) . This means that, very close 
to the diabolical point, the frequency l^k actually has a 
conical dependence on the wave number. Later we will 
understand the consequences of the existence of these 
points. 



The fact that 



1,2) 



: constant causes that 
do not have a wave character; only 

are true waves. Hence propagation in the 2D 
Grover walk occurs only if the initial condition projects 
onto subspaces 1 or 2; otherwise the walker remains lo- 
calized. This explains why strong localization in the 2D 
Grover map is usually observed, as first noted in Ref. 
liot - All this can be put more formally in terms of the 
group velocity Vg (k) of the waves in the system, given 
by the gradient (in k) of the wave frequency. The group 
velocity, as in any linear wave system, has the meaning 
of velocity at which an extended wavepacket, centered in 
k-space around some value ko, moves (there are other 
effects affecting extended wavepackets which we analyze 



below). In our case uj- 



(s=3,4) 



are constant, hence their 



gradient is null: these eigenvalues entail no motion. On 
the contrary cj^*"^'^^ depend on k (|27|28p and then define 
non null group velocities: 



,(1,2) 



(k) = iVkl^k, 



(30) 



4 



k2ln 0.0 




-0.5 



Figure 2: Vector plot of the group velocity (shown as arrows), 
superimposed to a density plot of the velocity modulus Eq. 

(EH}. 



where the plus and minus signs stand for s - 
spectively. Using (|28|) these velocities become 



,(1,2) 



(k)=± 



(sin fcijSin 



1^4 — (cos ki + cos ^2)^ 



1,2, re- 



(31) 



Hence the maximum group velocity modulus in the 
2D Grover walk is as can be easily checked and is 
obtained, in particular, along the diagonals /c2 = 
This is the maximum speed at which any feature of the 
QW can propagate. Figure 2 displays a vector plot of 
Vg ' (k) where a number of interesting points (in k 
space) is revealed: On the one hand there are points 
of null group velocity, located at k = (0, ±7r) and at 
k = (±7r,0), in correspondence with the saddle points 
observed in Fig. 1. On the other hand there are five sin- 
gular points at k = (0, 0) as well as at the four corners. 
These points are singular because the group velocity is 
undefined at them (they look like sources or sinks). We 
note that these points coincide with the diabolical points 
in Fig. 1. 

At this point it is worth mentioning that with coin 
operators other that the Grover diabolical points are also 
found, as we have checked by using the DFT coin [46| 
that displays several conical intersections. In this case of 
the DFT coin none of the four leaves of the dispersion 
relation is constant. We mention too that a behavior 
very similar to that of the 2D Grover walk is found in 
the Alternate QW j4o|, except for the fact that in this 
last case there are no constant surfaces in the dispersion 
relation (and hence no localization phenomena). 



A. Eigenvectors 

The eigenvectors of the 2D Grover walk matrix (jSj 
are given by Q 



4^) 





\- 








1- 








\ - 








'\ - 







(32) 



where 



n 



(5) 



-ifci 



)( 
)( 



1 



(33) 



is a normalization factor, and A^^^ = exp i^^^^^ 
is the associated eigenvalue. Use of these eigenvalues in 
([24j) allows finally to compute the state of the system at 
any time. 

Whenever k is not close to a diabolical point these 
eigenvectors vary smoothly around k. Here we just want 
to study the behavior of the eigenvectors close to the 
diabolical point at k = ko = (0, 0) -at the corners the 
conclusions are similar. We find it convenient to use polar 
coordinates (/ci,/c2) = (/ccos6>, /csin^). Performing the 
limit of (|32]) for /c ^ we find 



2V^ 



2V^ 



1 
71 



/ 1 + 

1 - 

-1 - 
-1 - 



1 - 
14 

-1 

-1 



V2cos(9 \ 
cos i9 

- V^sin(9 

- sin i9 / 

V^cos(9 
cos i9 

- V^sin^ 

- V^sin^ 





sin^ \ 








— sin6> 


6^^^ - - 


1 




cos6> 




1 


v 


— cos d 1 




VI/ 



(34) 



Hence, close to the diabolical point ko, there are three 
eigenvectors, corresponding to s = 1,2,3, displaying 
a strong azimuthal dependence, while (\)^ is constant 
around ko (its variations are smooth and hence tend to 
zero as k ^ ko). Thus, even if only ^^^^ and ^^^^ partic- 
ipate in the conical intersection and define the diabolical 



point. 



is also affected by this feature because it is 



associated with eigenvalue A 



(3) 



-1, which is degener- 



ate with A^*^^'^\ The fact that the diabolical point is 
singular is easy to understand: depending on the direc- 
tion we approach to it the eigenvectors of the problem 
are different. 



5 



it 



Just at the diabolical point A^^ ^'^'^^ = —1 and 
hence there is a 3D eigensubspace formed by all vec- 
tors orthogonal to the fourth eigenvector, namely 

^ col (1,1,1,1), which is equal to (j)^^: see ([34]) . A sensi- 
ble choice for these eigenvectors follows from noting that 
^ + ^k^^ in (|34|) is independent of the azimuth and is 
orthogonal to 0^"^^ Hence one of the basis elements for 
this 3D degenerate subspace is chosen as 



1 



V2 



col (1,1, -1,-1), (35) 



which has the outstanding property of being the single 
vector that projects, close to the diabolical point ko, 



only onto 



i.e. onto the eigenspaces of the dia- 



bolo. As for the other two eigenvectors we just impose 



their orthogonality with ( 
arbitrarily as 



and 



and we choose them 



■1,0,0), 



^ col (0,0,- 



■1,1). 



(36) 

Note that we are using a "primed" notation instead of 
keeping the notation with the label 5, except for 5 = 4, 
because none of these eigenvectors keeps being an eigen- 
vector for k ^ ko (only for special directions 6) and 
hence there is no sense in attaching any of these eigen- 
vectors to a particular sheet of the dispersion relation. 
To conclude we list the projections of these eigenvectors 
onto the eigenvectors (|34|) in the neighborhood of ko : 



1 

(2)\ cos6> 
k / - 72 ' 
(2)\ _ sin6> 



(3) 
k 



0, 



smt 



I ^L^^) 



and 



= cos^, 

(37) 
0. Ah this 



has strong consequences on the QW dynamics near ko 
as we will show below. 



IV. CONTINUOUS WAVE EQUATIONS FOR 
SPATIALLY EXTENDED INITIAL CONDITIONS 

In this Section we describe the evolution of spatially 
extended initial conditions that are close to the plane 
waves we have introduced. We define such spatially ex- 
tended states as those having a width which is apprecia- 
bly larger than the lattice spacing (taken as unity in this 
work). This type of initial states are wavepackets which 
are easily expressed in reciprocal space. 



2N 



(s) 

k-ko 



(38) 



where kg is a reference (carrier) wavevector, chosen at 
^k ^ / associated eigenvectors, and is 



will. 



a narrow function of k, centered at k = (i^^-ko ^ 
centered at k = ko) and having a very small width 
Ak^'^ < TT [aj. Notice that we have chosen an initial 
coin state which is independent of k. From here one 
must distinguish between regular points^ where eigenvec- 



is 



tors 



j have a smooth dependence on k close to ko, 
and diabolical points^ where eigenvectors have wild vari- 
ations around them, as we have seen. 



According to (p!3|) the initial condition 
position representation 



lV^x,o) 



2N 



=1 



reads in 



(39) 



where 



is) 



x,0 



/ 



{2n) 



N' 



,i(k-ko)-xj^(s) 



k-ko 



(40) 



~ (s) 

is a wide and smooth function of x because is con- 
centrated around k = 0. Hence in real space our ini- 
tial condition consists of a coin state (/>ko^^ equal at all 
points, multiplied by the carrier exp [i (ko • x)] and by a 
wide and smooth function of space, F^^q. We see then 
that the type of initial conditions we are dealing with 
are very close to the plane waves of momentum ko of the 



QW, given by e'^^'"^ 



A. Regular points 

For regular points the initial condition (|38]) determines 
the weight functions (|2T]) as 



k-ko 



O 



(41) 

As) 

"ko 



where we took into account that the eigenvectors 
vary smoothly around k = ko and that only for k 

~ (s) 

the function i^-ko nonvanishing. Note that we cannot 
make such an approximation in the case of diabolical 
points because of the strong variations of the eigenvectors 
around them. 



Now, the system will evolve according to (f24|) . Ap 



proximating the eigenvector 



appearmg m 



by 



using the same arguments as before, the partial 



6 



waves 



i^xl ) can be written as 



= exp 



(ko 



k-kn' 



where we have defined new functions -i F^^} | 



O (aA;(*)) , 
(42) 

(43) 

jFi'.U . Note 



^ exp (i [(k - ko) • X- (J"^ 

[271) ^ ^ 



that, at t = 0, V^^^Q 
agreement with (j39j). 



^i'o exp [i (ko • x) 



The problem is then solved if functions 



2Ar 



are 



s=l 



determined. Instead of doing so by brute force, i.e. by 
integrating -maybe numerically- Eq. (|43l) . what we do 
now is to look for a continuous wave equation that, by 
comparison with other known cases, sheds light onto the 
expected behavior of the QW. 

In order to avoid confusion, we introduce a function 
F^*^ (x, t) of continuous real arguments exactly in the 
same way as in (|43|) as there is nothing forbidding such 
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a definition. We have then F^^j. = F*^^^ (x, t) for x G 
and t G N. In order to simplify the derivation we rewrite 
Eq. (|43|) by making the variable change k — ko ^ k. 
Finally, as for the limits of integration (originally from 
— TT to +7r for each dimension), we extend them from — oo 
to +00 in agreement with the continuous limit. We then 
have 

nv exp 



/ N f r / 

(x,t) = / -—^exp 1 (k-x 



t 



XF 



(44) 



Let us obtain the wave equation. First we take the 
time derivative of (|44|) and get 

d^k 



idtF^'^ (x,t) 



(2^) 

■^ko+k 



AT 



exp 



i k • X- 



^ko+k - ^ko 



t 



(45) 

Then we Taylor expand the function I — ^ko (ex- 

cept in the exponent: otherwise large errors at long times 
would be introduced) around k = 0, 



N 

^ko+k-^ko - Z^^i 

i=l 



1 



N 

2! ^ 



where 

is) _ 



k=0 



k=ko 



^ = d'uj^^ljdhdkj 



k=0 



dJ^^/dk, 
- d'^J^^ /dkidh 



(46) 



(ko) 
(47) 



k=ko 



(48) 



etc. and we have recalled that Vg^^ (ko) = Vk^^ 

Now it is trivial to transform the so-obtained right hand 
side of Eq. (|^ as a sum of spatial derivatives of 
F^^^ (x, t), leading to the result 



k=ko 



(x,t) 



dt 



(ko)-VF(^)-| E 



dxjdxn 



(49) 

which is the sought continuous wave equation. One 
should understand that, in general, few terms are neces- 
sary on the right hand side of Eq. (|49]) because F^^^ (x, t) 
is a slowly varying function of space, as commented. 

The first term on the right hand side of the wave 
equation is an advection term, implying that the ini- 



tial condition Fq^^ (x) = F^^^ (x, 0) is shifted in time as 
F(^) (x, t) = F^'^ (^x - v^^^ (ko) without distortion, to 
the leading order. In fact we can get rid of this term by 
defining a moving reference frame X such that 

X = x-v(^)(ko)t, ^(^)(X,t) = F(^)(x,t), (50) 

and then Eq. ()l9|) becomes 



' dt 



1 f „<..^ 

2 dXidXj 

1,3 = 1 



(51) 



which is an A/'-dimensional Schrodinger-like equation, to 
the leading order. This means that the evolution of 
the wavepacket consists of the advection of the initial 
wavepacket at the corresponding group velocity and, on 
top of that, the wavefunction itself evolves according 
to (|5T]) . Equation (|5T]) is a main result of this article 
as it governs the nontrivial dynamics of the QW when 
spatially extended initial conditions, as we have defined 
them, are considered. It evidences the role played by 
the dispersion relations as anticipated: For distributions 
whose DFT is centered around some ko, the local varia- 
tions of uj around ko determine the type of wave equation 
controlling the QW dynamics. 

Just two additional remarks: First, if the initial condi- 
tion only projects onto one of the sheets of the dispersion 
relation (hence the sum in s in (|38|) reduces to a single 
element), ah F^^) (X,t) wih be zero but the one corre- 
sponding to the chosen eigenvector. This means that 
the coin state is preserved along the evolution, to the 
leading order. Second, if the initial state projects onto 
several eigenvectors the probability of finding the walker 
at (x, t), irrespectively of the coin state, follows from (j5j), 
(|23l) and (02]), and reads 



p.,=^|Fi:) =^U(^)(x,t) 



(52) 



to the leading order, where we used (^^^j | ^) = 5s^s'- 
Hence Px,t is just the sum of partial probabilities: there 
is not interference among the different sub-QW's. 
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1. Application to the 2D Grover walk 



In order to illustrate how the evolution of the QW 
can be controlled via an appropriate choice of the initial 
conditions, we will show some results obtained numer- 
ically from the evolution of the 2D quantum walk with 
the Grover operator Eq. (|25]) . Unless otherwise specified, 
the initial profile in position space is a Gaussian centered 
at the origin with cylindrical symmetry (independent of 
a or T] ): 

V^x,o=Are^^°-"e-4^, (53) 

with X = (xi,X2), which implies that in the continu- 
ous limit we also have a Gaussian, in momentum space, 
centered at ko = (/ci,/c2)- Here, A/" is a constant that 
warranties the normalization of the state. In the numer- 
ics, we have taken a sufficiently large value <j = 50 for 
the Gaussian, so as to make it possible the connection 
with the continuous limit, although we will also discuss 
the situation for smaller values of a in order to show the 
robustness of the results. 

Let us first consider a simple case having ki = 
k2 = 7r/2, with an initial coin state ^(| (/>ko^)+ | 

cPI^^)), where | = -1= col (1, 0, -1, 0) and | 0^)) = 

^ col (0, 1, 0, —1) that projects on the 5 = 1 and s = 2 
branches with equal weights. The group velocity is given 
by Vg^^ (ko) = ±(^, ^), respectively for s = 1,2, imply- 
ing that the original Gaussian will split into two pieces 
that will move along the diagonal of the lattice in oppo- 
site directions, as shown in Fig. [3l It must be noticed 
that, since the second derivatives vanish at this point 
{ujij = 0), the Gaussian distribution moves with no ap- 
preciable distortion. On the right bottom panel of this 
figure, we show the result for <j = 5 and we observe that, 
in this case, the results obtained for large a are very ro- 
bust even for such a low value of a. 

However, if we choose an initial value kg close to the 
origin, then the second derivatives take a large value, 
with the consequence that now the two pieces experience 
a considerable distortion along the line perpendicular to 
the line of motion, as can be seen from Fig. |4l where 
^1 = ^2 = O.OItt, implying a curvature vuij o:^ 7.96 in the 
perpendicular direction. 

In this case, at variance with the case considered in 
Fig. [31 we have to face a peculiar situation, as the distri- 
bution is now peaked around a value close to the diabol- 
ical point ki) = (0,0). As far as the initial distribution 
does not contain this point with a significant probability 
(i.e., for a large value of a), the behavior will be as de- 
scribed above. However, if a decreases below a certain 
limit, the evolution is dominated by the singular behav- 
ior of the diabolic point (to be studied in detail in the 
next subsection). This transition is clearly observed in 



-4 
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Figure 3: (Top panel): Probability distribution as a function 
of the dimensionless (xi,X2) position. The initial condition 
is given by ^(| 4>^^^)+ \ and fei = fca = 7r/2.. The 

distribution is shown for three times t=0, 400, 800 respec- 
tively. (Bottom panel): Top view of the previous figure, both 
for cr = 50 (left) and a = 5 (right). The numbers indicate, in 
each case, the value of time for different snapshots. 

Fig. [5l where the resulting evolution is depicted as a is 
increased from cr = 5 to cr = 30, the latter subpanel al- 
ready showing a close resemblance with the result in Fig. 
m Distributions initially centered around non-singular 
points, however, are much more robust against a smaller 
cr, as previously shown in Fig. [3l 

A completely different situation arises when ko cor- 
responds to one of the saddle points (such as the one 
located at k^p = (0,7r)). As the group velocity van- 
ishes at that point, i.e., ' (ko) = (0,0), the cen- 
ter of the probability distribution is expected to remain 
at rest but will spread in time analogously to diffract- 
ing optical waves, as Eq. (|5T]) reduces to the paraxial 
wave equation of optical diffraction [72]. In our case, 
by considering 5 = 1 as the initial coin state, one finds 
wii = — l/2,tz7i2 = tU2i = 0,tz722 = 1/2. Therefore, Eq. 
(f5T]) becomes 

. aA(^)(x,t) _ l aM(i) 

' dt "4 dXf ^4 dXl ' ^ ' 

This wave equation is similar to that of paraxial opti- 
cal diffraction [72], differing from it in that the sign of 
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Figure 5: Same as Fig. [H corresponding to four different 
values of a. From left to right and from top to bottom: a = 5, 
cr = 10, cr = 20 and a — 30, respectively. The initial condition 
is given by M\ 4>\^h+ I C)) ^1=^2 = 0.01^. 




Figure 4: (Top panel): Probability distribution as a function 
of the dimensionless (xi, X2) position. The initial condition is 
given by ^(| )+ | )) and /ci = /cs = O.OItt. The dis- 
tribution is shown for three times t=0, 400, 800 respectively. 
(Bottom panel): Top view of the previous figure, showing the 
value of the corresponding time steps. 



spatial derivatives are different for the two spatial direc- 
tions. Hence this hyperbolic equation describes a sit- 
uation in which there is diffraction in one direction and 
antidiffraction in the other direction, as it happens in cer- 
tain optical systems [73] . This causes that the cylindrical 
symmetry proper of diffraction is replaced by a rectan- 
gular symmetry. Basing on the analogy with diffraction, 
it is possible to tailor the initial distribution in order 
to reach a desired asymptotic distribution as we already 
demonstrated in the one-dimensional QW [l3|. Hence 
if we want to reach a final homogeneous (top hat) dis- 
tribution, we must use an initial condition of the form 



'sinc{x/(Jo) sinc{y /(Jo) ^ (55) 



with sinc{x) = sin (ttx) /ttx and ctq a constant that ac- 
counts for the width of the distribution. In Fig. [6] 
we observe the evolution of the Grover walk with the 



above initial condition and the initial coin 4^^^ 0^ = 

^co/ (1 + 1 + i, 1 — i, 1 — z), equal for all populated 
sites. The figure shows a final distribution that is quite 
homogeneous along most of its support, except in its 
outer borders, where it shows a smooth but rapid fall 
out (the ratio a/ao has been chosen to optimize the result 
0). This is the expected result, but Fig. [6] also shows 
a central peak that has been cut for a better display, 
the height of this central peak being quite large (around 
1.5 • 10~^ ) as compared to the plateau (around 4 • 10~^). 
The existence of that central peak is a consequence of 
the common initial coin |^ko,o)- This coin state guaran- 
tees that the initial distribution projects just onto the 
relevant branch of the dispersion relation at its central 
spatial frequency but, as the width of the distribution is 
finite, for the larger values of |k — ko| it will provide a 
non-negligible projection onto the static branches of the 
dispersion relation and the localized part of these projec- 
tions are what we see as a central peak in Fig. [6l 

B. Diabolical points 

This case requires special care because there are eigen- 
vectors of the QW that vary strongly around the con- 
ical intersection, forbidding the very initial approxima- 
tion taken in the case of regular points, namely Eq. (|4T]) . 
Given the singular nature of the problem we will try to 
give a rather general (but not fully general) theory here: 
We will present a treatment that covers the case of the 
2D Grover walk and, very likely, can be applied to other 
cases. The basic assumptions are: (i) There are just 
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1000 



Figure 6: Probability distribution as a function of the di- 
mensionless (xi,X2) position at time t=9000, for initial con- 
ditions at the saddle point ko = (0,7r). The coin initial state 
is I 0k '*) — I col (1 + i, 1 + i, 1 — i, 1 — i), and the starting 
space distribution as given by Eq. ([55)) . We have taken a 
value (Jo = 15 for this simulation. The central peak has been 
cut for a better display of the rest of the features. 



two (hyper-)sheets in the dispersion relation that dis- 
play a diabolical point at k = ko, and we label them 
by 5 = 1,2 for the sake of definiteness; and (ii) There 
are other sheets, degenerate with the diabolical ones at 

(s) 

k = ko, whose associated frequencies uj-^^ are constant 
around the diabolical point. 

In order to alleviate the notation we will denote by uj^^ 
the value of the degenerate frequencies at k = ko- In 
particular 



(1) (2) 



but there can be other sheets with 

,(3) 



,(^7^1,2) 



(56) 
uj^ (as it 

happens with ijj)^^ in the 2D Grover map). 

We consider an initial wavepacket defined in reciprocal 
space by 



'k,0 



Fx 



k-kr 



|5) 



where |S) is an eigenstate of the coin operator Ck at 
k = ko, and Fk-ko is again a narrow function centered 
at k = ko- In the position representation this initial 
condition reads, see (p!3|) . 



V^x,o) 



,0 IS) 



(58) 



where Fx,o = / 



d^^k ^ifk-k 



^^""^k-kB is the DFT of F^^. 



Hence, as in the regular case, |S) represents the state of 
the coin at the initial condition. Clearly |S) is not unique 
because of the degeneracy. We do not use the notation 
^Ld) instead |S) because the former are ill defined, 
as we have seen in the example of the 2D Grover walk. 



The initial condition (|57|) determines the weight func- 
tions ([2T]) as 



-As) 



i^k- 



(59) 



We will assume that |S) does not project on regular 
sheets (those not degenerate with the conical intersec- 
tion); otherwise those projections will evolve as described 
in the previous case of regular points. This is equivalent 
to saying that, in the remainder of this subsection, all 
sums on s will be restricted to sheets that are degenerate 
at the conical intersection. 

Definition (|59|) poses no problem but for k = ko . How- 
ever, as this is a single discrete point (a set of null mea- 
sure) its influence on the final result, given by integrals 
in k, is null [74] . 

A main difference with the regular case is seen at this 
stage because there is not any choice of the initial coin 
state |S) making that only one value of s be populated 

(remind that ^k'') ^^^^ strongly around k = ko). This 
has consequences as we see next. 

Now the system will evolve according to (|24|) . The 



partial waves ) read now 



gi(kD-x-u;j3t) 
d^k 



I 



(27r) 



N 



exp 



I 



kx- 



6^^) ) 



(60) 



(61) 



(57) E 



where we did the variable change k — ko ^ k, and ex- 
tended the limits of integration to infinity as correspond- 
ing to the continuous limit. Note that ^k^^'^^ — ujj^ = 
according to assumption (ii) above. As in the regular 
case, the evolution has a fast part, given by the carrier 
wave e^(^^"^~'^D^) ^ and a slow part (both in time and in 
space) given by F^^^^^ 

Note that the previous approximations reproduce the 
correct result at t = 0, Eq. (|58|) , upon using that 



'^kD+k/ \V-kD+k 



is very approximately the unit op- 
erator in the degenerate subspace, to which |S) belongs. 

It is important to understand that the vector (coin) 
part of the state evolves in time in this diabolical point 
situation, because more than one of the sheets that be- 
come degenerate at the conical intersection become nec- 
essarily populated, i.e. there are at least two s values for 



^ko+k 



S) 7^ 0, unlike the regular case. 



which 

Finding a wave equation in this case is by far more 
complicated than in the regular case, as the vector part 



of 



depends on time (because 



cannot be 



approximated by 
it was the constant 



kD+k/ 

), unlike the regular case where 
Then we will not 
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try deriving such a wave equation but conform ourselves 
with trying to understand the evolution of the system 
under the assumed conditions. 

What we can say is that the frequency offsets 



00- 



(s=l,2) 
kD+k 



(f6T]) are conical for small k (i.e. in 
the vicinity of the diabolical point, which is the region 



,(5=1,2) 



""kD+k 



selected by Fk); in other words, 

where c is the group speed (the modulus of the group 
velocity) and k = |k|; see the 2D Grover walk case in Eq. 
(|29|) . where c = 1/a/2- For the rest of sheets we have 
assumed cj^^^^'^^ — = 0- Hence we can write, from 

m 



1. The 2D Grover walk 

In this case N = 2, d^Q = d6>, with 6 the azimuth, 
and 

E(?^ (fc)) = d^exp (ik . X) (^1;^^^ I S) I ^1;^^^). 

(66) 

We remind that here |S) is assumed to project just onto 
^ki^k^)- "^^^^ means that \E) coincides with |0d), see 
(35|). According to the projections (|37j) and to the vectors 
(|34|) we can write, in matrix form. 



— — ^ exp [i (k • X T ckt)] E^^^ (^) = 7 / d(9 exp (ik • x) 

(27r) 4 Jo 



I -) I ^£+k), 



(62) 



while the rest of partial waves verify 
/-^exp[i(k.x) 



r(s/i,2) 



^kD+k/ 



x,0 



(63) 



which is a constant, equal to its initial value. Hence, the 
projection of the initial state onto the sheets of constant 
frequency does not evolve, as expected, leading to a pos- 
sible localization of a part of the initial wave packet. We 
thus consider in the following that (^kj+k ^ I ^) ~ ^• 

Regarding the sheets 5 = 1, 2 we cannot progress unless 
we assume some property of the eigenvectors ^^D+ky 
As in the 2D Grover walk case, we will assume that 
^ki^k^) (for k close to zero, i.e. around the diaboli- 
cal point) depend just on the angular part of k (let us 
call it Q) but not on its modulus k. This does not seem 
to be too a strong assumption. We also restrict ourselves 
to cases in which Fk only depends on k (spherical sym- 
metry), i.e. F\^ = Fk. Hence we write the integral in (|62]) 
in (TV-dimensional) spherical coordinates as 

Fi7'''^^ = i^^y^ J dfcexp(TicH) k^-^F^ \e[^'^ (fe)) , 

(64) 



E(f) (k)) = I d^l^exp (ik • x) {cPI:^^^ I -) I 0^^^^^,), 



(65) 



where we wrote d^k = k^~^dkd^Q and d^Q is the A^- 
dimensional solid angle element. 

Up to here the theory is somehow general. However we 
cannot progress unless we particularize to some special 
case. Next we do this for the 2D Grover walk. 



/ l±V2cosl9 \ 
1 T V2 cos l9 
-1 T V2 sin l9 
V -l±V2smO ) 



(67) 

where the upper (lower) sign in ± or ^ corresponds 
to s = 1(2), respectively, and we remind that k = 
{kcosO^ks'mO). Upon writing k-x = kxcos{6 — (p), 
where x = x (cos (/p, sin cp) and performing the integral 
we get 



E(,«) (fc) 



/ Jo (kx) ± i>/2 Ji (kx) cos (p \ 
Jo (kx) ^ iV2Ji {kx) cos cp 
— Jo (kx) ^ iV^ Ji {kx) sin (p 
\ — Jo {kx) ± i>/2 Ji {kx) sin (p ) 



(68) 



where Jn are the Bessel functions of the first kind of order 
n. Using this result in (|64|) we get, in matrix form, 

F^'] = (27r)"^ ^ dA:exp {^ikct) kFj^E^^^ {k) . (69) 

Let us compute the probability of finding the walker 
at some point, regardless the coin state, as given by (j5j), 
([23]) , (|6Q1) and (|69|) . After simple algebra we get 



it 



where 



Px^t — (2'^) ^ / d/ccos {kct) Jo {kx) kPk^ 
Jo 

POO 

p^^l = (27r)"^ / dk sin {kct) Ji {kx) kFk, 
Jo 



(70) 

(71) 
(72) 



are amplitude probabilities. We see that the spatial de- 
pendence of the probability amplitude on the azimuth (p 
in ()69l68p disappears when looking at the full probabil- 
ity Px,t- This means that a coin-sensitive measurement 
of the probability should display angular features. 

As a summary of what we know up to here about the 
neighborhood of the diabolical point, we can say that 
when the initial coin |S) coincides with |0d), see ([35]) . 



11 



we predict (i) no localization and, (ii) when the initial 
wavepacket Fx,o is radially symmetric (just depends on 
X = |x|), an evolution of the probability Px,t that keeps 
the radial symmetry along time, according to ()7Q|71|72p . 

Equations (fTT]) and ([72]) describe the evolution of the 
probability of an initial wavepacket centered at a diabol- 
ical point and with radial symmetry. The actual proba- 
bility Px,t can be computed from them numerically but 
little can be said in general. In order to gain some insight 
we consider next the relevant long time limit, which turns 
out to be analytically tractable. First of all it is conve- 
nient to scale the radial wavenumber k to the (loosely de- 
fined) width of the initial state in real space, Fx,0 5 which 
we denote by a. Accordingly we define = ak so that 
(fTT]) and ([72]) become 



We readily obtain 



(27r(j^) ^ J dA^cos (~^) ^"^^^ 

!} = {27Ta') Ak sin (^^k^ Ji (^«) (74) 



where = F^^^/^ is non null only ioi n < 1. Hence 
when ct ^ cr^ cos (^^) and sin (^^) are strongly oscil- 
lating functions of a^, around zero, what would make p^^'P 
to vanish. However in the integrals defining such am- 
plitude probabilities other oscillating functions, Jo (^/^) 
and J I (^Ai:), appear. Thus it can be expected that when 
the oscillations of the latter and the ones of the former are 
partially in phase, a non null value of the integrals is got. 
According to the asymptotic behavior of (2:) at large 
Jn(z) ^ Y^2/7r2: cos (z - n7r/2 - 7r/4) for z > - 1/4| 
[75|, this partial phase matching occurs when x ^ ct, 
which is expected from physical considerations. Hence 
we are led to consider the limit x = ct + a^, with ct^ (j^ 
where ^ is the scaled radial offset with respect to ct, in 
which case ([73|) and ([74|) . after expressing the products 
of trigonometric functions as sums, become 



(0) (1) 



(27r)"Vv'2^ / d/csin(/ccr^ + 7r/4) V^Ffc, (75) 
^0 



^(0) _ (1) _ 



1 



(76) 



where ^ = and In is the modified Bessel func- 

tion of the first kind of order n. The total probability 
Px,^ = P(0 follows from using (|76|) in ([70]) . from which 
we observe that the initial width cr acts only as a scale 
factor for the height and shape of the probability at long 
times (ct ^ a). A plot of the probability can be seen in 
Fig. [71 where two unequal spikes, whose maxima are lo- 
cated at ^ = —1.74623 and at ^ = 0.550855, separated by 
a zero at ^ = —0.765951, are apparent. We remind that 
X = ct + cr^ is a radial coordinate, hence these maxima 
correspond to two concentric rings separated by a dark 
ring. The situation is fully analogous to the so-called 
"Poggendorff rings" appearing in the paraxial propaga- 
tion of a beam incident along an optic axis of a biaxial 
crystal, a situation in which a conical intersection (a dia- 
bolical point) is present. Remarkably the result we have 
obtained for Px,t fully coincides with that described in 
(59I whose Fig. 7 is our Fig. [71 

As an example that illustrates the above features, we 
have studied the evolution of a state of the form Eq. ([53]) 
for the diabolic point ko = (0, 0), and a coin state as de- 
fined in Eq. ([35]) . As seen from Fig. [S] the probability 
distribution keeps its cylindrical symmetry while it ex- 
pands in position. It can also be seen from this figure 
that the general features remain valid for a wider distri- 
bution (in k space) corresponding to a = 5, even if the 
distribution in spatial coordinates is narrower for a given 
time step. The Poggendorff rings can be appreciated in 
Fig. [9] which shows a radial cut of the previous figure. 

All these features crucially depend on the choice of the 
coin initial conditions, as clearly observed when these 
conditions are chosen differently. As an example, we 
show in Fig. [10] the evolution after 200 time steps, of 
the probability distribution when we start from a state, 
also centered around the diabolic point, but the vector 
coin is now ^^d) ^^^^ 6> = 7r/2. Most of the probability 
is directed along the positive X2 axis, which corresponds 
to the choice of ^. In fact, by changing the value of this 
parameter one can, at will, obtain a chosen direction for 
propagation. 



to the leading order, where the approximation follows 
from disregarding highly oscillating terms. As an appli- 
cation of this result we consider next the Gaussian case 
Fk = 2cry^exp (— ^a^/c^), where a is the standard de- 
viation of the initial probability in real space, Eq. ([53]) . 



V. CONCLUSIONS 

In this article we have formally solved the standard 
multidimensional QW taking into account extended ini- 
tial conditions of arbitrary shape and width. In the limit 
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Figure 7: The probability P(^), scaled with the prefactor act, 
as obtained form Eqs. (|7Q|) and (|76|). 




500 -500 




Figure 8: (Top panel): Probability distribution as a function 
of the dimensionless (xi,X2) position, at time t = 400. The 
initial condition is given by Eq. ()35|) around the diabolic point 
ko = (0, 0). (Bottom panel): Top view of the previous figure 
with a = 50 (left) and a = 5 (right), but now at time t = 100. 

of large width the continuous limit is particularly well 
suited and clearly reveals the wave essence of the walk. 
The use of the dispersion relation is the central argument 
behind this view as its analysis allows to get much insight 
into the propagation properties of the walk, even allow- 
ing for the tailoring of the initial distribution in order 
to reach a desired asymptotic distribution, as we have 
demonstrated. We have also shown that the two dimen- 




Figure 10: Probability distribution as a function of the di- 
mensionless (xi,X2) position, at time t = 200. The initial 
condition is given by l^k^) ^tround the diabolic point with 
6* = 7r/2 . 

sional Grover QW exhibits diabolical points in its disper- 
sion relations and have analyzed in detail the dynamics 
around this point, closely connecting the Grover walk 
with the phenomenon of conical refraction (Hsl, [60| . Al- 
though we shall not give details here, we have also found 
conical intersections when the DFT coin is used instead 
of the Grover one, but in the DFT case there are more 
conical intersections and the structure of the dispersion 
relation (that lacks of constant surfaces) is quite more 
involved than in the Grover case. 

Our results have also shown that the asymptotic dis- 
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tributions found in the continuous limit, valid for large 
initial distributions, can be qualitatively reached for not 
very broad distributions, say a 5, and this should be 
easy to implement in systems such as the optical inter- 



ferometers of Ref. [28|, |37| . An exception to this general 



rule are those distributions which are peaked close to the 
diabolic point, since in this case a broad distribution will 
be dominated by the singular nature of that point. 
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